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Flow in both saturated and non-saturated vuggy 
porous media, i.e. soil, is inherently multiscale. 
The complex microporous structure of the soil 
aggregates and the wider vugs provides a multitude 
of flow pathways and has received significant 
attention from the X-ray computed tomography 
(CT) community with a constant drive to image at 
higher resolution. Using multiscale homogenization, 
we derive averaged equations to study the effects of 
the microscale structure on the macroscopic flow. The 
averaged model captures the underlying geometry 
through a series of cell problems and is verified 
through direct comparison to numerical simulations 
of the full structure. These methods offer significant 
reductions in computation time and allow us to 
perform three-dimensional calculations with complex 
geometries on a desktop PC. The results show that the 
surface roughness of the aggregate has a significantly 
greater effect on the flow than the microstructure 
within the aggregate. Hence, this is the region in 
which the resolution of X-ray CT for image-based 
modelling has the greatest impact. 



1. Introduction 

The macroscopic Darcy's law and Richard's equation 
which describe flow in porous media can be derived 
either using formal two-stage homogenization [1-4], 
two-scale convergence [4,5] or volume-averaging techni- 
ques [6,7]. Homogenization has been widely used to 
describe flow in single porosity materials [1,4,8,9] and, 
more recently, to describe flow and diffusion in dual 
porosity models where the substructure is composed 
of two different domains with different porosity 
[10,11]. In these structures, Darcy's law is applied in 
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each domain and homogenization allows an averaged Darcy law for single-phase flow or 
Richard's equation for two-phase flow to be derived [10,12,13]. 

Another related structure is vuggy porous media [14-17]. Vuggy porous materials consist of 
regions of tightly packed microparticles which form aggregates separated by larger pores or vugs. 
In the vuggular region, Darcy's law is no longer applicable and Stokes flow must be considered 
[15]. The macroscopic behaviour of flow in such a medium relates not only to the flow in the vugs 
and the aggregate, but also to the condition applied at their interface. This condition is geometry- 
dependent and is often described using the Beavers and Joseph condition [15], or the Saffman 
approximation to the Beavers and Joseph condition [18]. Both these conditions are slip boundary 
conditions which relate the shear stress to the velocity at the surface of a porous medium with 
slip length proportional to the square route of the hydraulic conductivity in the porous medium, 
where the constant of proportionality is left as an experimentally determined fitting parameter. 

There have been several studies which aim to eliminate the fitting parameters used in the 
Beavers and Joseph conditions [14,19-22]. Levy & Sanchez-Palencia [14] derived simplified 
boundary conditions for two limiting cases: the case where pressure gradients were normal 
and non-normal to the interface. The rigorous derivation of the Beavers and Joseph condition, 
which generalizes these cases, was first derived by Jager & Mikelic [19]. This method uses 
the assumption that the porous domain is periodic in the direction tangential to the boundary 
enabling the Beavers and Joseph coefficient to be derived through a cell problem relating the 
average velocity in the Stokes domain to a unit shear stress applied at the boundary. 

Recent studies by Arbogast et al. consider vuggy porous geometries [16,17]. Arbogast et al. 
apply the Saffman approximation of the Beavers and Joseph condition on the boundary between 
the vug and the adjacent porous medium [15,18]. This condition is based on the fitting parameter 
in the Beavers and Joseph condition which depends on the exact geometry of the boundary 
considered. The result of these studies is a macroscopic derivation of Darcy's law in which the 
hydraulic conductivity depends on the coupled flow in the vugs and the aggregate. 

In this paper, we extend the work of Arbogast to include the geometrical properties of the 
interface between the vugular region and the porous region. Specifically, we study the flow 
of fluid in vuggy porous media in the context of two-phase flow in soils. In order to answer 
fundamental questions regarding flow in porous media and the interaction of these flows with 
external sources and sinks, e.g. roots, it is essential to develop a model which captures all 
necessary geometrical features of the soil [23]. Not only will this model provide significant 
insight into the flow mechanisms and advanced models which can be incorporated into image- 
based simulations [24], it will also feedback into the resolution driven imaging of soils through 
X-ray computed tomography (CT) and synchrotron radiation-based microtomography [25,26] by 
providing a lower limit to the scale of soil features which affect flow properties and hence, need 
to be detected by X-ray CT. 

We consider the flow of air and water in a periodic array of soil aggregates (figure la). The 
aggregates are composed of a periodic array of soil particles (figure Ih). The resulting geometry 
has three different scales: the soil particle scale, the aggregate scale and the macroscopic or field 
scale. This structure has been designed to have a bimodal pore size distribution as observed 
in measurements on typical soils [27-29]. The hydraulic properties of this geometry are highly 
complex. To simplify the hydraulic properties of the fluids, we consider the case where the 
aggregate is completely hydrophilic and the interface between the two fluids is stationary [4]. The 
dynamics of fluids in vuggy porous media are highly complex [30-33]. These geometries exhibit 
circulation near interface between the vug and the porous region even in the case of low Reynolds 
numbers [33], this is attributed to flow penetration from the vug into the porous region [34,35]. 
However, these effects are seen to decay over time and settle to the steady-state flow properties 
[33]. We consider the commonly used steady-state Stokes equations, [4,16,17], and neglect the 
short time scale dynamics of the fluid. To understand the precise role of geometry in such a 
structure, we apply the boundary conditions of Jager & Mikelic [19-22]. Our resulting model 
is free from fitting parameters and describes the flow in vuggy porous media defined entirely by 
the geometry and fluid properties. 
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Figure 1. Idealized soil schematic showing typical length scales and pore sizes, (a) Shows the macroscopic picture, a periodic 
arrangement of aggregates and air bubbles, (b) shows a zoomed in image of the aggregate scale, this is a single unit cell of 
the macroscale geometry showing the internal structure of the aggregate and (c) shows a zoomed in image of the microscale 
geometry inside the aggregate. (Online version in colour.) 

The application of this theory requires the assumption of periodicity of the pore structure on 
the aggregate surface [19], something which is not easy to achieve even in idealized geometries. 
However, the results for hydraulic conductivity are shown to be accurate in the cases tested. 
Typical errors for the geometry shown in figure 1 are 2% for the Beavers and Joseph condition 
and 10% for the Saffman condition for a highly porous aggregate with pore sizes ^23% the 
size of the maximum space between aggregates. The errors reduce to 1% for the Beavers and 
Joseph condition and 2% for the Saffman condition in the case of a low porosity aggregate 
with typical pore size of ~ 8% the size of the space between the aggregates. In non-ideal soil 
geometries, the assumption of periodicity is less likely to be accurate. However, this can be 
addressed by considering successively larger sample sizes and imposing periodicity. As the 
sample size increases the errors induced by imposing periodicity are likely to reduce and the 
macroscopic properties will converge. The resulting hydraulic conductivity depends on the 
underlying periodic structure of the aggregate, the nature of the flow at the aggregate surface and 
the flow around the aggregate itself. It is shown that the required boundary conditions depend on 
the relative scales of the inter-aggregate pores, the vug radius and the macroscopic sample size. 

This paper is arranged as follows: in §2, the theory is derived and for flow in vuggy porous 
media. In §3, the boundary conditions tested for two- and three-dimensional idealized aggregates. 
Finally, the results are discussed with reference to imaging in §4. 



2. Model 

The derivation of the model is arranged as follows: first, in §2a, we derive the scaled equations 
governing fluid flow in the porous medium. In §2^, we apply homogenization theory to the 
equations and study the behaviour of flow in an individual soil aggregate. In §2c, we study the 
boundary conditions at the aggregate surface in terms of the theory developed by Jager & Mikelic 
[19,20]. Finally, in §2d, we derive the macroscopic equations which describe the behaviour of the 
averaged flow in the soil while taking the microstructure into account. 



(a) Scaling 

We consider the flow of two incompressible fluids contained within the porous medium (figure 1) 
with length scales and fluid properties applicable to the flow of water and air in soil [36,37]. 
The soil consists of a periodic collection of aggregates (figure la) with an internal microstructure 
which is also periodic (figure lb). We refer to the flow around a collection of aggregates as the 
aggregate scale and the flow through the pore structure within the aggregate as the microscale. 
The structure is porous on both the aggregate scale and the microscale with typical pore size 



Ty and fz, respectively. The aggregate itself (figure lb) is assumed to be roughly periodic both 
with distance around the surface and internally This assumption is clearly not completely 
valid even for simple geometries such as this one and for simulations carried out on X-ray CT 
geometries it is necessary to increase the size of the periodic unit considered until the hydraulic 
properties converge. The structure has fundamental period Ly on the aggregate scale and on 
the microscale. The macroscopic length scale, Lj (figure la) is chosen based on the macroscopic 
geometry considered. Here, we consider the case of Lj ^ 10 cm, a length comparable with the 
typical length scale for soil columns used in X-ray CT imaging [26]. The maximum inter-aggregate 
pore size is typically ry^lOOi^m [29]. In unsaturated soils with no air sources or sinks, the 
dominant force is gravity. For gravity-driven Poiseulle flow, we find that the Reynolds numbers 
for the water and air phases are 

Rei-^J^P4<1 and Re^^^ = ( (2.1) 

where /x^^^ and /x^^^ are the water and air viscosity, p^^^ and p^^^ the water and air density and g 
is the acceleration because of gravity. Here, we concentrate on highly tortuous porous media and 
as such we expect this to be an over-estimate of the Re3molds number. Hence, we consider the 
Stokes limit of the Navier-Stokes equations. 

We also assume that the system is in dynamic equilibrium such that the location of the 
boundary between the two different fluid phases is known/ fixed. We denote the water domain 
^2w/ the air domain ^2^, the soil particle boundary Tg and the air-water interface Taw The 
dimensional Stokes equations, in an unsaturated geometry [4], are 

X e Q^, (2.2a) 
X e Q^, (2.2b) 

xeQysj (2.2c) 

X e ^2a, ('2.-2d) 

where a^^) = (Vu^"^^) + (Vw^'^^)^ and a^^) = (Vu^^^) + (Vw^^^)^ are the stress tensors for the fluid 
phases, p^^^ and p^^^ are the water and air pressures, respectively, u^^^ and ii^^^ are the water 
and air velocities, S^^^ and S^^^ are the water and air saturation which satisfy S^^^ + S^^^ = 1 and 
0 is the porosity denoting the total fraction of volume available for flow. We note that equations 
(2.2c,d) are macroscopic equations which refer to the saturation of the soil. The rigorous derivation 
of these equations requires careful consideration of the fluid-fluid interface dynamics and is 
beyond the scope of this paper. The macroscopic equations we have used are commonly accepted 
in the homogenization of two fluids in porous media [4] and, as we will show in §2^, on the 
microscale reduce to the standard incompressibility conditions for each fluid. The porosity, 0, is 
also a macroscopic quantity and refers to the total fraction of space available for flow of either 
fluid. Assuming periodicity of the structure and that the soil matrix itself remains stationary then 
0 does not change over space or time. On the soil particle surface we use a no slip condition 
combined with zero fluid penetration. Hence, all the velocity components vanish on the surface: 

zi(^)=0, M^^^^O, xeTs. (2.2e) 

We also define a set of boundary conditions on the air-water interface. Specifically, we require 
that the interface is stationary, i.e. the normal velocity of the two phases is zero 

u^^^'h = 0, w(^^-w = 0, xeTaw, (2.2/) 

the slip length associated with tangential stress goes to zero for Stokes flow, hence, the tangential 
velocity is continuous 

[^(w)_^(a)].^.^Q ^^^^^^ ^2.2g) 



and 



95(w) 
dt 



+ V . u^^^ = 0, 



and there is a jump in the normal stress given by the surface tension curvature product 

n . [/x^^^a^^^ - J^(^^] - h . Ifi^^^a^^^ - ip^^^] = hky xe Taw (2.2/z) 

Here, k and y are the curvature and surface tension of the air-water interface, respectively; h and 
f y, for ; = {1, 2}, are the vectors normal and tangent to the air-water interface, respectively, with h 
pointing into the water domain. We scale space to the fundamental period on the aggregate scale 
(figure lb), x = Ly^x such that the aggregate scale is periodic with period 1. We define two small 
parameters, e = Ly/Lx denotes the ratio of the aggregate scale to the macroscopic length scale 
and T] = Lz/Ly denotes the ratio of the microscopic length scale to the aggregate length scale such 
that, locally to the aggregate, the microscale is periodic with period r]. We introduce the following 
non-dimensional variables: 



and the dimensionless parameters are 

which we use to derive the non-dimensional equations prior to using multiple scale 
homogenization to derive the equations macroscopic which describe the flow of air and water 
in soil. This scaling results in equations which are effectively the same as the ones used in [4]. 
The key differences are that we have chosen to scale the viscosity and density of the two fluids 
into the parameters 8p and 8u and that we have chosen to scale with the aggregate length scale 
Ly. This choice of spatial scaling results in a slightly different expansion of the gradient operator, 
we consider variations on the scale small rather than considering variations on the Ly scale 
large. However, this choice seems a more natural as it makes it easier to keep track of the different 
spatial scales. The non-dimensional Stokes equations which results from this scaling are 

e V . a(^) - V]y(^) = eg, xe Q^, {23a) 

e V . a(^) - Vjy(^) = e^g, X e Q^, {23b) 

a5(w) 

60 + V . M^^^ = 0, X e (2.3c) 

dt 

a 5(a) 

and 60 + 5w V • u^^^ = 0, xe ^2a, {23d) 

dt 

with boundary conditions 

m(^) = 0, u^^^ = 0, xeTs, {23e) 

u^^^'n = 0, M^^^-w^O, xeTaw, (2.3/) 

W^]-fy = 0 xeTaw {23g) 

and h • - Ip^"^^] - 8ph . [ea^^^ - Ip^^^] = nK xe Taw (2.3/z) 

Here, the non-dimensional stress tensors are given by cr^^^ = (Vm^^^) + (Vm^^^)^, cr^^'> = {Vu^^^) + 
(Vm(^))T. 

We note that the scaling of the air velocity and pressure is different to the scaling on the 
water velocity and pressure by a factor of e. This is justified both physically and mathematically. 
Physically, it can be seen using equation (2.3^) that it takes an air velocity of order to induce a 
water velocity of order 1. Mathematically, it can be seen that scaling the velocity of air and water 
identically would lead to the replacement of 8u by €8u in equation (2.3^) and 8p by €8p in equation 



(2.3/z). The result of this change is that it is impossible to balance the equations at 0(1) when the 
homogenization procedure is applied. 

The e in front of the scaled gravitation term in equation {23h) comes from the difference in 
scaling of the water and air velocities. This choice of scaling is appropriate for the flow of air and 
water in soil. However, the method described here is widely applicable for Stokes flow and it is 
trivial to adapt the resulting equations for different gravitational scaling. 

(b) Microscale behaviour 

We start by considering how the microscale geometry affects flow on the aggregate scale. In order 
to do this, we seek an asymptotic solution to the scaled equations (2.3). We do this first in terms 
of e to isolate the flow behaviour in the aggregate scale unit cell (figure lb), then in terms of 
T] to isolate the flow behaviour in the microscale unit cell (figure Ic). Initially, we define two 
different spatial scales x denotes the macroscale spatial coordinate, y = ex denotes the aggregate 
scale spatial coordinate. Expanding the velocity and pressure in powers of e we obtain 

= u^^\x, y) + €U^^\x, y) + 0(£^) and j^^"^ = p^^\x, y) + €p^^\x, y) + 0{e^), (2.4) 

where a = {w, a} denotes the water and air phase, respectively, and apply the two-scale 
homogenization procedure [3], V = Vy + £ V^. Collecting terms at order £^ we obtain 

Vyp^^^ =0, ye Q^, {2.5a) 

'^yPo^ = 0. yeQ^ {23b) 

and p^^^ - 8pp^Q^ =k, ye Taw (2.5c) 

fa) (a) 

Equations (2.5) have solution, pr^ ^jr^ {x), i.e. the leading order pressure is a function of the 
macroscopic variable only. Furthermore, equation (2.5c) relates the macroscopic air and water 
pressures through the surface tension curvature product. In order to determine the effect of 
the macroscopic pressure gradients on the microscale flow, we proceed to the next order of 
the expansion and obtain a set of equations for the leading order velocity and the order e 
pressure correction 





y e ^2w, 


{2.6a) 




y e Qy^, 


{2.6b) 




yeQ^, 


(2.6c) 




yeQ^, 


{2.6d) 






(2.6c) 


(w) ^ i^i) n 

u-Uq =0, fi'uy = 0, 


y e Aw, 


(2.6/) 


Tj '{uy - Suu'q 0 = 0/ 


y e Aw, 


(2.6g) 


n-l{at^-Ip^r^)-8,{4'^-Ipt^^^ 


y^r^w 


(2.6/2) 


p^^^ and Wq"^ are periodic with period 1, 




(2.60 



where = {Vu^) + {Vu^)^ and = {Vu^ + {Vu^^. At this point, we could solve 
equations (2.6) to obtain the macroscopic hydraulic conductivity of the soil. However, as the 
aggregate is made up of small particles with tiny pores, this would be computationally intensive. 
To overcome this, we search for approximate solutions to equations (2.6) which average out the 
flow in the aggregate. 

As the structure in the aggregate is periodic with period r], we introduce the new variable 
z = y/ri to denote the microscopic spatial coordinate and expand the gradient operator as 



Vy = {l/ri)Vz + Vy. As a result, we have three different spatial coordinates capturing the three 
different spatial scales x, y and z to denote spatial position on the macroscale, the aggregate scale 
and the microscale, respectively. 

For simplicity, we assume that the aggregate is completely saturated with water. Hence, we 
need only consider one-phase flow in this region. We denote the aggregate domain Q^, the extra- 
aggregate water domain as ^2ws = -^w \ expand the Darcy domain velocity and pressure 
in rj 

^cT^ = ^ao(^' V' ^) + ^^04 V' ^) + ^^^fl^^' V' ^) + V^^d (2.7fl) 

and 

v'x^ = Vx^{^.v,^) + nvf^{^,v,^) + n^vfx^^,y.^) + 0{n\ yea^, (2.7b) :| 

with the fluid equations 

vy-^ - Vyp^-^ = Vj-^ ye (2.8.) 

Vy . u^^^ = 0, ye Q^, {2M) 
t/[r^ = 0, yeFs (2.8c) 

and and Mq^^ are periodic with period 1. (2.8d) ' 

We are interested in writing an averaged equation which describes the behaviour of the fluid 
on the microscale. The microscale geometry may be considered periodic assuming that we are 
far from the aggregate boundary. To obtain an approximation for the flow in this domain, we 
substitute equations (2.7) into (2.8) and consider a single periodic unit cell on the microscale 
(figure Ic). Collecting the terms in equation (2.8) in ascending powers of 7], we obtain Wq^q = 0, p^^^ 
is invariant on the microscale, i.e. p^^^ = p^^Q{x,y) and u^^^ = 0. We note that Mq*^q = 0 and u^^^ = 0 
effectively rescales the Darcy velocity to 0(r]^). However, in order that the velocities internal and 
external to the aggregate remain scaled the same we retain the original scaling. We expand to 
order ri^ to obtain the well-known result for single-phase Darcy flow, see for example [4], with an 
additional source term owing to the x dependence of the pressure gradients 



3 

^0^^ = E n{z)[dyj^^{y) + d,J^\x) -^h'g] (2.9a) 



and 

3 

PU = E 0^k{z)[dyj^^{y) + d,J^\x) + h ' g]. (2.%) 

k=l 

Here, is a unit vector in the kth direction and and oj^ satisfy the cell problem 

V^vfc - V.ojk = h, ye (2.10a) 

Vz • Vfc = 0, ye Q^, (2.101?) 

vfc = 0, yeTs (2.10c) 

and co'^Q and Vq^-^^ are periodic with period r]. {2.\{)d) 

Expanding equation (2.8^) to 0{r]), integrating over the Darcy domain, applying the divergence 
theorem and using equation (2.9), we derive the averaged equation for the microstructure 

^yK^ypfl = Q, (2.11) 

where K is the hydraulic conductivity of the aggregate, 

Ci • vy dy. (2.12) 



X2d 



This result allows us to apply Darcy's law on the interior of the aggregate with the microscale flow 
driven by the macroscopic pressure difference. However, in order to obtain the total macroscopic 
flow in the medium, we need to consider flow on the aggregate scale. Before we can do this, 
we need to define appropriate boundary conditions for the flow at the interface between the 
aggregate, in which Darcy flow holds, and the adjacent pore space, in which Stokes flow holds. 



(c) Boundary condition 

In order to complete the approximation for the hydraulic conductivity on the aggregate scale, 
we consider the boundary condition on the aggregate surface. We follow the method of Jager & 
Mikelic [19-22]. 

The key observation is that the Darcy velocity in the aggregate scales with r^, where Vz is the 
pore radius while in the pore space the velocity scales with r^, with Vy the radius of the macropores 
external to the aggregate. Typically, we expect that the ratio of these numbers r^/fy ^ L^/Ly = r] 
such that u^^^ ^ 0{r]^) and m^^^ ^ Using this assumption, we define the rescaled average 
aggregate velocity 

ri^i,^^ = -<Pd'K{Vyp^^l + V,p<"') + g), (2.13) 

where 0^ = is the microscale porosity defined as the total fraction of the representative 

microscale volume available for flow. In the previous section, we derived Darcy's law for flow 
through the aggregate. This allows us to write the equations on the aggregate scale as 









(2.14fl) 


Vy 




y e Qws, 


(2.Ub) 


Vy 






(2.14c) 


Vy 






(2.14d) 


Vy 




yeQd, 


(2.14e) 



and 

where the boundary conditions (2.6/-/z) still hold. We now determine the behaviour of the flow at 
the interface between the aggregate and the adjacent pore space. We assume that the aggregate 
surface is strongly hydrophillic and, as such, the air-water interface remains at a distance from the 
interface which may be considered large on the microscale. This allows us to derive the boundary 
conditions assuming that only water is in contact with the aggregate. We define a false interface 
denoted S and, as we did in §2^, expand the fluid velocity and pressure outside the aggregate in 
powers of r] 

^0^^ = ^ao(^' ^) + ^^04 2// ^) + ^^Wa2(^/ ^) + y e ^2ws (2.15a) 



and 



p^f = pflix, 1/, z) + r^pflix, t/, z) + rj^pflix, i/, z) + 0{r^^l y e Q^s- {2.15h) 



Our aim is now to match the Darcy velocities and pressures to the pore space velocities and 
pressures of equal order in r] on S. Substituting (2.15) into (2.14) and collecting terms of order rf^ , 
we obtain the Stokes problem for the velocity and pressure of the water and air phases 

Vy'^o!^^-Vypi"o^-V,p[,")=0, yeQ^s, (2.16fl) 

Vy«ao^=0, yeQ^s, (2.16fe) 

Vy^ao-Vyp'^^-V4^^ = 0, yeQ, (2.16c) 

and Vyu[,^^ = 0, yei2a, (2.16d) 



where a'^' = (Vyu''^^) + (V,u<'^^)'^ aj^^ = (VyU^^l) + (V^u^^lf 

(w) 
^0,0 





y ^ J, 


(2.16e) 




y e Taw, 


(2.16/) 






(2.16^) 


0, 


ye^aw 


{1.16h) 



(w) ^ - (a) n 

n • Wq Q = 0, n • Mq Q = 0, 

/ (w) o (a) X ri 

and « . [(a<;^> - Ip^^o^) - - 

The assumption of no slip at the interface induces a jump discontinuity in the shear stress which 
is equal to h - a^^^ip, for p = {1, 2}. The presence of a shear stress at such a false interface is clearly 
unphysical. To correct for this, we consider the behaviour of the water in the region close to 
the aggregate surface. Specifically, we consider a boundary layer of width ^ 0{ri) and, as in §2^, 
rescale y = r]z, Vy = V^, and look for a velocity and pressure which balance the shear stress 
jump in this region. Expanding the velocity and pressure in the boundary layer we find that the 
velocity balances at order r] and the pressure at order 1: 

M^l = ^M^l + _^ Q(^3) ^bl^^bl_^^^bl_^Q(^2)_ (2.17) 

The boundary layer problem which results from expansion (2.17) is 







z e Qy^, 


(2.18fl) 






z e ^2w/ 


(2.18^7) 




«^' = o, 




(2.18c) 




«^'ls+-«^'ls-=0. 


zeS, 


(l.lSd) 






zeS 


(2.18e) 


and 


P^, periodic with period 77 in ti, f2direction 




(2.18/) 


with (Tq ' = (VzMq') + (VzMq')^. Equations (2.18) are separable; hence, we can write the velocity and 
pressure local to the interface as 




2 2 


(w) / X ^ 


(2.19) 


where ;8^' 


and Xp^ satisfy the cell problem 








v,a^'-v,x^' = o. 


z e ^2w/ 


(2.20fl) 






z e ^2w/ 


(2.20^7) 






zeFs, 


(2.20c) 




<ls+-<ls-=o. 


zeS, 


{220d) 




n • [(a^l - Ix^%^ - (a^l - Ix^%-] = -Tp, 


zeS 


(2.20c) 


and 


^ , Xp^ periodic with period ^ in ti, f2 direction 




(2.20/) 



and ^ = (^zPp^) + (^zPp^)^- Here, S+ denotes the water side and S~ denotes the aggregate side 
of the false boundary and n is a normal vector pointing out of the aggregate. As equation (2.20) 
is defined on an infinite domain in the direction normal to S, the solution method is non-trivial. 
Results from [21] conveys that the velocity tends to zero with distance from the false boundary in 
^rid that the normal velocity tends to zero with distance from the boundary in However, 
the tangential velocity tends to a constant in Furthermore, the pressure has non-zero average 



at the boundary which will induce an additional flow in the aggregate. We deflne the far fleld 
velocity, C^J, and the pressure jump, C^J^, as 



pbl 



^T,-p^'dz and Cl = 



xfdz. (2.21) 



s 



In order to solve these equations on a flnite-sized geometry, Jager et al. [21] propose a set of 
additional constraints which enforce the predicted far field behaviour of the flow. Specifically, 
as the velocity tends to zero with distance from the false boundary in we apply a no-slip 
condition at a finite distance, L^, from the false interface. In the Stokes domain, the normal velocity 
tends to zero with distance from the false boundary and the tangential velocity tends to a constant. 
Hence, we apply a slip boundary condition at a distance Lw from the false boundary 

P'f = 0, z e Ld, {1.11a) 

h-p^^ = 0, zeL^ (1.11b) 

and h-G^^Tq = 0, with = {1,2}, zeLw (2.22c) 

As these conditions enforce the predicted behaviour, we expect that as and Lw increase the 
solution will converge. Typically, this happens for L^ = L^>lr]. 

As jSp^ decays to a non-zero constant with distance from S, the solution extends into the 
bulk and must add an additional contribution to the flow. This counter flow is generated in 
the aggregate scale pore space along with an additional pressure term at the boundary which 
contributes to the flow in the aggregate. Specifically, the additional flow satisfies Stokes equations 
in the water and air domains 

• ''ti - ^yvti (2.23fl) 
Vy«[,^> = 0, yeQ^s, (2.23b) 

• '^ai - '^yPu f ^ (2.23c) 

and Vyu^^=0, yei2s, (2.23d) 

where a^J^ = i^yU-^i) + (^yWo^V and a^f = (Vyu[,^J) + (^yU-^})^. The additional velocity at the 
false boundary must be matched to the velocity in the pore space 

«0A «al=«o'' l/eTd (2.23e) 

and the first-order correction to the air-water boundary conditions in r] are 

(w) ^ - (a) n 

n ' Uq-^ =0, n ' Uq [=0, 

/ (w) o (a)>. p. 
Tj • (M^y -5i,M^j) = 0, 

and ^.[(,(j)_j^(w))_,^(,(^^^^ 

In the aggregate Darcy's law, derived in §2^, holds with the boundary pressure given by the 
pressure in the fluid domain and the additional correction from the boundary layer 





y e Taw, 


(2.23/) 






(2.23g) 


0, 




(2.23/z) 



and 



Ptl = Vt^ + Vti + E ^l^ • -t^'v y ^ ^d, (2.24^.) 



where we recall that the microscale hydraulic conductivity K is defined by equation (2.12). The 
next order correction comes from the Darcy velocity at the boundary. As the Darcy velocity is 





{2.25a) 




(2.25b) 




(2.25c) 


zeS 


(2.25d) 


zeS, 


(2.25e) 



small, KVyp^^Q ^ 0{r]^), it will induce a correction ^ 0{r]^) in the boundary layer. We write the 
correction to the boundary layer problem as follows: 

= 0, 

,M\ „bl| „('^) 
Ml \S+ - 1*1 Is- = 1*0,2' 

and h ■ [(fTi" - Ip\%. - (ai" - Ip1%-] = h ■ V.u'^^ 

where a^^ = (V^m^^) + (V^m^^)^. As in the leading order boundary layer problem we could solve 
equations (2.25) to obtain the velocity near the boundary and consider the counter flow generated 
in the pore space. However, as the modulated part of the solution decays with distance from the 
boundary and the constant part is zero in the Darcy region we expect that the velocity in the pore 
space will stabilize to the average velocity at the interface of the porous medium. 

The average velocity at the interface is simply the Darcy velocity divided by the aggregate 
porosity 0^/ see equation (2.13). Therefore, the 0{r]^) counter flow will be generated directly by 
the Darcy velocity. We note that the stress contribution on the particle scale is the integrated 
average of the normal derivative of the tangential velocity on the microscale. This is identically 
zero and the stress jump from the flow in the Darcy region does not directly contribute to the flow 
in the pore space. The counter flow generated by the Darcy velocity is given by 
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and the boundary conditions 




(w) -(d) (a) -(d) 
^0,2 ~ ^av / ^0,2 ~ ^av / 






1/e^d, 


(2.26c) 




(w) ^ ^ (a) n 
n • Mq 2 — 0/ W • Mq 2 = 0, 






y e Aw, 


(2.26/) 




^; • (>2 ~ ^w>2) = 0' 
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{2.26g) 


and 


^ r/ (w) T (w)x e> / (a) 




= 0, 


yeFaw 


(2.26/2) 



Rather than solve the hierarchy of cell problems which generate the flow in the aggregate 
and the adjacent pore space, it is computationally more efficient to use the behaviour in the 
boundary layer to write an averaged boundary condition on the aggregate scale. By combining 
the contributions at different orders, we can write the boundary conditions on the false boundary 
as follows: 

h . u^^^ = ri^h . uf^ {2.27a) 

2 

fp • u^^^ = n J2 • ^r^^) + ri\ • u^^ {227b) 



and - ) = e ^ Cl{h • a^^f,). (2.27c) 

Clearly, these conditions are only correct to 0{r]^) and more terms must be considered if higher 
accuracy is needed. 



As we are only interested in the leading order solution to the homogenized problem, the 
maximum order of terms considered in the boundary condition must depend on the ratio -q/e. 

li r)^ 0{e), then the boundary layer terms correspond to a correction order e and need not 
be considered at this order. If this is the case then the appropriate condition on the aggregate 
boundary is the no-slip condition. For the case where rj^ ~ 0{e) then the first-order correction in r] 
must be considered and the appropriate boundary condition is the Saffman condition. Finally, if 
T]^ ^ 0{e), then the appropriate boundary condition is the original Beavers and Joseph condition. 
We will revisit this point in §4. 



(d) Aggregate scale 

We now return to the equations for the macroscale hydraulic conductivity, i.e. the flow on the 

X scale (figure \a). Using the aggregate scale hydraulic conductivity from §2^ and the boundary 
conditions from §2c, we can write equations (2.6), for all orders in rj, as 

' 4"^ - ^yvT^ = VxFo"^ +Sr ^ws, (2.28a) 

Vy . = 0, ye ^2ws, (2.28^7) 

. ^(5^^ - \y^f = y e ^2a, (2.28c) 

Vy . u^^^ =0, yeQ^ {228d) 

and VykVypf^ = 0, yeQ^, (2.28e) 

with the boundary conditions at the air-water interface 

h-u^^'>=0, h'U^^^ = Q, yer^y,, (2.28/) 

ij • {u^^^ - 8uu^^^) = 0, ye Taw (2.28g) 

and h . [{ea^''^ - Ip^^^) - S^iea^^^ - Ipf^)] = nK, ye Taw, (2.28/z) 

and the Stokes-Darcy interface 

h • u^^^ = r^^h . k{^ypf + ^,p^^^ +gl ye Tdw, (2.28/) 

2 

. = ^ ^ C)\{h . a^^^ . f^) + . U'^ypf + ^,p^^^ +gl ye T^w (2.28;) 



2 

and pf^ - p^^^ = e ^ C^l (n • a^"^^ • f,), y e T^w, (2.28^) 

P=i 

where rj^K = (p^^K and K ^ 0(1). Equations (2.28) are separable and, hence, the solutions can be 
written as 

3 

= E '^^"''"\y){dxj^\x) + h-g) + ic^^''\y)d,J^\xl with a = {w, a} (2.29) 



and 



3 

Pi^ = E '"^ W(9x.F[r^ W + h-g) + 4''''\y)d,J^\x), with a = {w, a, d}, (2.30) 



where tt^'^ and k^'^ satisfy the cell problems for ^ = {a, w} 

Vy4'^'^> = 0, i/ei2ws, (2.31fc) 

Vy$f'^>-Vy7r*"'^> = Sa/;efc, i/ei2a, (2.31c) 

Vy • = 0, ye Qa, (2.31d) 

Vy • kVyjtl'^'^^ =0, y e i2d (2.31e) 

and n^"''^^ and /c*"'^^ are periodic with period 1, (2.31/) 

with = (Vy,.["''^)) + (Vy/.["''«)T, 5^^'^) = (Vy/.[^'^)) + (Vy/.[^'«)T and the air-water boundary 
conditions 

H.4"''«=0, «.4"'^) = 0, ye Taw, (2.31^) 

rp.(4"''^>-S„4^'^)) = 0, i/eTaw (2.31/1) 

and ^i.[(5(w,«_j^(w,«)_g^(-(a,«_^^(a,«)^^0^ y^^^^^ (2.310 

and aggregate boundary conditions 

h • K^^'^^ = rj^h • k(Vy7Tjf'^^ + S^^hl y e Tdw, (2.31;) 

2 

Tp ' k^^^ - ^'^(VyTrf + 8^ph)] = r) J2 • 4"^'^^ ' y ^ ^dw (2.31^) 

and nf'^^ - tt^^^ = ^l^^ ' ^^''^ ' V ^ ^dw, (2.31/) 

where Soifi = 1 for a = ^ and ^q;^^ = 0 otherwise. The additional source terms in equations (2.31;,/c) 
come from the homogenization procedure on the aggregate surface. Physically, these terms 
provide the hydraulic conductivity contribution owing to the macroscopic pressure gradient 
across the aggregate. We also obtain through the first-order correction and application of the 
Fredholm alternative the macroscopic Richard's equation for saturated flow 

r\ C(W) 

<l>^ = V,(/C(-'-)v4"" + K;("'^)v4^') (2.32a) 

and 

cl>^= S„V,(;c(^'-)V,p|,"') + >C(^'^)V,p|f)), (2.32fe) 

where icf/^ = - 4^'« • Cj dy and ^^'^ = - J^(4-'« + .^KVy^rf '«) • Cj dy. 

Equations (2.32), the capillary pressure equation (2.5c) and the saturation condition S^^^ + 
5(a) _ 2 combined with the cell problems described in §2h,c describe the pressure-driven 
saturation of soil in a vuggy porous medium with the hydraulic conductivity parameters 
determined entirely based on the aggregate geometry. 

This model captures the flow of two fluids in and around the aggregate for the case in which 
the aggregate is strongly hydrophilic. This model can be simplified in the case of single-phase 
flow. In this case, S^^^ = 1 and we can see that equations (2.31) simplify to give a single cell 



problem for a:^^) and 



(w,w) 



V~ (w,w) (w,w) /V ^ ^ /n oo \ 

Vy . XVyTif '"^^ = 0, yeQ^ (2.33c) 

and TT^^'^^ TT^^'^) and /f ^^'^^ are periodic with period 1, (2.33^i) 

with the aggregate boundary conditions 

n • K^^'"^^ = n^h • kiVyjzjf''^^ + h). y e r^w, (2.33.) 

f , . [^-'-^ - ^^X(V,;rf + a,)] = , C)\{h . . f 1/ E (2.33/) 

q=l 

2 

J (d,w) (w,w) V^/^bl/^ ~(w,w) ^ s. /o oo \ 

and H ^^k =1^ * ^fc * '^v)' y ^ ^dw (2.33^) 

The result of this simplification is that /C^^'^^ = /C^^'^^ = /C^^'^^ = 0 and, neglecting terms 0{rl^), 
equation (2.32) reduces to the standard Darcy's law for vuggy porous media [10]. However, we 
should note that the constant in the Saffman approximation is derived from the geometry rather 
than left as a fitting parameter. Similarly, in the limit that the porosity of the aggregate goes to 
zero, we can recover the two-phase flow models described in [4] simply by neglecting terms 0{r]) 
in equations (2.31). 

The theory derived in this section highlights the role of the different space scales in the 
derivation of the macroscale hydraulic conductivity in soils. Physically, the relative sizes of r) 
and £ conveys the macroscopic length scale at which the assumptions made on the aggregate 
boundary are applicable within the approximations made in the aggregate scale homogenization 
procedure. The macroscopic length scale can be written L^f^ = Ly{Ly/LzY , where L^^^ is the length 
scale on which the no-slip condition may be used on the aggregate surface and the higher order 
terms do not contribute significantly to the flow. The length scale tells us the length at which 

Saffman's simplification of the Beavers and Joseph boundary condition produces a contribution 

(3) 

to the flow of order 1, L\ ' is the length scale at which the fully coupled model must be solved. 
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oo 



3. Example 

In this section, we numerically study the effect of the boundary condition used on the macroscopic 
hydraulic conductivity. We consider two cases both of which are axially symmetric about the 
'false' boundary normal. In this case it can be proved that the pressure correction C^J^ = 0 [21]. In 
§3fl, we study the hydraulic conductivity of a two-dimensional system with different aggregate 
structures. In §3^, we consider the effect of the aggregate structure on a three-dimensional 
idealized soil sample. 



(a) Two-dimensional geometry 

We consider two different two-dimensional geometries using the method illustrated in figure 1. 
The first geometry consists of a periodic packing of circular aggregates with radius 0.35Ly. 
The second geometry we have considered for the two-dimensional case is a square aggregate 
(figure 7.g). The square is aligned to the principal coordinate axis, has side length 0.4Ly and 
has smoothed corners with radius of curvature 0.05Ly. The microscopic aggregate structure is 
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Figure 2. Numerical solutions showing absolute dimensionless velocity, [a) Shows the microscale cell problem, (b) the 
numerical solution with the geometrically derived Beavers and Joseph boundary condition, (c) the full numerical solution with 
the aggregate structure taken into account, (d) shows the solution for an aggregate with a solid core and (e) shows the solution 
for an aggregate with a hollow core, (f) shows the solution in the boundary layer with velocity streamlines and (g) shows the 
full numerical solution for the square aggregate geometry. (Online version in colour.) 



also periodic and is composed of particles which are ellipses with principle radii Vj =RxLyrj 
and Tn =RnLyr] and {RmRj} e [0,0.5] with Rj =0.5 or Rn = 0.5 corresponding to the particles 
touching. Here, the particle scale is ^ = 0.05, the resulting minimum micropore radius is given 
by min{(l - lRn)r)Ly/l, (1 - lRr)r]Ly/l} . 

The cell problem for aggregate hydraulic conductivity, (2.31), depends only on the parameters 
S-p and 8u • These are determined by the viscosity and density of the two fluids and the ratio of 
the macroscopic- to-microscopic length scales. In our simulations, we have used /x^^^ = lO^^ Pa s, 
/x(^) = 1.8 X 10"^ Pa s, yo(^) = 10^ and p^^^ = 1.2. The typical aggregate size is Ly = 10"^ m [37], and 
we consider a macroscopic length scale of Lj = 10~^ m, chosen to be comparable to length scales 
for soil columns used in X-ray CT imaging [36]. The resulting dimensionless parameters are bp = 
8.47 and 8u = 0.15. 

There is clearly a question as to the applicability of the boundary condition on the curved 
surface of the aggregate and, for the circular aggregate, the periodicity of the particles inside. 
Here, we assume that the boundary conditions derived in §2c are appropriate and that the errors 
induced by this assumption will be small. 

We calculate the boundary conditions and the macroscopic hydraulic conductivity for a variety 
of different-sized and -shaped particles within the aggregate using Comsol Multiphysics. The 
Saffman problem decouples and the Stokes problem is solved in isolation before the Darcy 
velocity is calculated in the aggregate. For the Beavers and Joseph case, the equations are solved 
iteratively in r/. First, the Stokes problem is solved, then the output is used as a boundary condition 
for the Darcy problem, this is then used to calculate the correction to the Stokes problem. We do 
this for four cases: the case where the full geometry is taken into account, the case where a no-slip 
condition is applied to the surface of the aggregate, the case where only the Saffman condition is 
applied and the case where the full Beavers and Joseph condition is applied (figure 2). In addition, 
to verify the role played by the internal aggregate structure, we have calculated the hydraulic 
conductivity without approximation for both aggregate shapes with a solid core and with a 
hollow core (see figure 2d,e for the geometry in the circular case). The full geometry case involves 
solving equations (2.6), the various approximations involve solving equations (2.31) neglecting 
terms of order rj for the no-slip case, r]^ for the Saffman case and order rj^ for the Beavers and 
Joseph case. 
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Figure 3. Numerically calculated macroscale dimensionless hydraulic conductivity for the full model, the case with the 
Saffman approximation and the case with the full Beavers and Joseph approximation in (a) two-dimensional circular, (b) two- 
dimensional square and (c) three-dimensional circular cases. (Online version in colour.) 

Table 1. Calculated values for C^' in the two-dimensional case for elliptic particles with principle radii Rj and 
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The boundary slip tensor C^J has size 1 is given for a variety of different ellipsoidal particles 
in table 1. It is clear from these results that the radius in the direction normal to the boundary has 
significantly more effect than the radius tangential to the boundary Physically, this is expected as 
the boundary layer problem is calculating the response of the tangential velocity to the surface 
geometry Increasing the particle radius in the normal direction reduces the size of the flow 
pathways within the aggregate resulting in a large decrease in the flow. Conversely, increasing 
the particle radius in the tangential direction leaves the critical flow pathway unchanged and has 
relatively little effect on the flow. 

The results for the macroscale dimensionless hydraulic conductivity are plotted for the square 
and circular aggregates in figure 3 for the circular inter-aggregate particles. These results show 
that the successive approximations for the boundary conditions produce higher accuracy results. 
As the microparticle radius is increased and, hence, the pore size is decreased the results converge 
towards the no-slip value. In the limit that the micropore size tends to zero, there is zero 
connectivity through the aggregate and the Darcy velocity must vanish, note this is not the case 
in three dimensions. The result is that the Saffman and the Beavers and Joseph conditions will 
converge to the same solution. 

Comparison between the aggregates with different internal geometries (figure 2c-e) shows 
that the internal structure of the aggregate has little effect on the hydraulic conductivity for larger 
particle sizes (figure 3) this trend is observed for both the square and circular aggregates. This is to 



be expected based on the accuracy of the Beavers and Joseph and the Saffman approximations at 
larger particle radius where the internal structure of the aggregate only comes into play at 0{7]^). 

The no-slip approximation which neglects all internal structure of the aggregate provides a 
poor approximation of the fluid properties. This improves with decreasing pore radius; however, 
it still performs poorly in comparison to the Beavers and Joseph or the Saffman conditions. In 
the case of the circular and spherical aggregates, the Saffman approximation, which includes 
the effects of the 0{ri) boundary layer, is significantly better than the no slip and reproduces 
the trends in the hydraulic behaviour as the microparticle size changes. It also provides a good 
quantitative approximation for the hydraulic conductivity. The Saffman approximation provides 
a less accurate approximation in the case of the square aggregate which is particularly notable for 
particle radius less than 0.25. This is because of the sharp corners present in the square geometry 
which induces a high-pressure gradient in the Stokes region. This results in a large Darcy velocity 
on the aggregate corner which decays rapidly with distance into the aggregate but induces a 
sizable counter flow in the Stokes domain at the corner of the aggregate. This counter flow is 
neglected in the Saffman approximation, as it comes into play at 0(77^). However, it is clearly 
needed to capture the effects of the corner at high aggregate porosity. 

The fully coupled case in which the flow in the aggregate is computed using Darcy's law 
coupled to the flow in the external pore space provides the best approximation. The boundary 
conditions and behaviour inside the aggregate are derived assuming periodicity of the structures. 
While the microscale inside the aggregate is clearly periodic the assumption of periodicity on the 
boundary is less accurate. The approximations used in the derivation require a flat interface with 
a periodic structure of microparticles. The interface curvature and the conflicting requirements 
for both a periodic structure inside the aggregate and on the aggregate surface leads to errors in 
the approximation. However, these are small even in the worst case, ^ 2% for the Beavers and 
Joseph condition. 

There is a clear advantage of using approximate boundary conditions rather than studying 
the full geometry in terms of computational speed. The full Beavers and Joseph approximation to 
the boundary condition is also clearly more accurate than the Saffman condition for small particle 
radius. However, the consequence is that the fully coupled model has to be solved for flow around 
and within the aggregate. 

(b) Three-dimensional geometry 

The advantage of approximation techniques is clearly greater in three-dimensional geometries 
where computation times can be large. Here, we consider the three-dimensional extension of the 
aggregate considered in §3a. The aggregate is spherical with radius 0.35Ly, the microparticles 
are spherical particles of radius RLyt] and again 7] is taken as 0.05. The internal structure of the 
aggregate is formed from spherical particles packed to form a cubic lattice with a shell, two layers 
thick, of spheres around the outside. Packing the particles evenly within the shell proves to be 
a much harder problem as there is no possible particle arrangement which provides an even 
distribution of spheres. However, in order that these methods provide a reasonable description of 
fluid flow in non-ideal geometries, it is essential that they are insensitive to small changes in the 
periodicity at the surface. Therefore, we have made no attempt to minimize the inhomogeneity 
in the sphere distribution and show numerically that the resulting hydraulic conductivity is 
still accurate. 

The macroscopic hydraulic conductivity is calculated for a range of different microparticle 
sizes, for each of the different cases considered in the two-dimensional example, using Comsol 
Multiphysics for the approximate cases and OpenFoam for the case in which the full geometry is 
taken into account. The equations implemented in OpenFoam are solved using the SIMPLE (Semi 
Implicit Method for Pressure-Linked Equations) algorithm [38], with an added source term. As 
the particles are spherical the diagonal components of the stabilization tensor are equal, = 
the off diagonal components are zero C\\ = C^J = 0. The value of C\\ is given for different R in 
table 2. 
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Figure 4. [a) Typical solutions to three-dimensional boundary layer problem and [b) aggregate scale cell problem with the 
Beavers and Joseph boundary condition showing absolute dimensionless velocity and velocity streamlines. (Online version 
in colour.) 



Table 2. Calculated values for C\\ in three-dimensional case for spherical soil particles of radius /?. 
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We now proceed as in the two-dimensional case and calculate the hydraulic conductivity 
for the four different cases. Typical solutions to the cell problems are shown in figure 4 and 
the results of these calculations are shown in figure 3. Both the Beavers and Joseph and the 
Saffman conditions offer a significant improvement in the estimated hydraulic conductivity when 
compared with the no-slip condition and highlight the importance of the approximate modelling 
techniques used in this paper. In terms of computational time, the approximate solutions took 
approximately 5 min to calculate for each radius on a desktop PC. The calculation which takes the 
full geometry of the aggregate took 20 h for each radius value running on 32 nodes of the IRIDIS 
high-performance computing facility at the University of Southampton. 

For low particle radii, corresponding to large pore sizes, the three-dimensional approximation 
behaves poorly with ^ 20% error for R = 0.1. This is not unexpected as, owing to the large pore 
sizes, there is significant flow through the aggregate and the errors induced by the disordered 
sphere packing and the surface curvature will become significant. However, it can be seen that 
for large particle radii, corresponding to small pore size, the error is significantly reduced to ~ 9% 
at R = 0.15 and ^ 2% for R = 0.35. 

The accuracy of the approximation techniques for large radius values, corresponding to small 
pore sizes, tells us that the internal aggregate structures is largely irrelevant in determining the 
macroscopic hydraulic conductivity. Rather, it is the pore structure on the surface of the aggregate 
that makes the most difference to the macroscopic properties of the structure. Hence, it is the 
surface roughness rather than the microstructure inside the aggregate which should be the focus 
of imaging techniques in order to obtain accurate solutions from image-based models. 



4. Summary 

In this paper, we have used the method of Jager & Mikelic [19-22] to derive the Beavers and 
Joseph boundary condition applicable to the surface of a soil aggregate in a periodic geometry. 
The resulting equations show how the hydraulic conductivity properties on a macroscopic scale 
relate to the geometry on the particle scale and the aggregate scale. 

The results show that the surface roughness of the aggregate is the key property of the 
microscale geometry which determines the hydraulic conductivity of the macroscopic geometry. 



Hence, this is the region in which the resolution of X-ray CT for image-based modelling has the 
greatest impact. More accurate approximations which take into account the coupling between 
the flow in the aggregate and the flow in the extra-aggregate pore space produce a slight 
improvement in the results at the expense of an increase in computation time, typically the 
Beavers and Joseph simulations take twice as long as the Saffman approximation using the 
iterative scheme mentioned in §3. 

By considering the different scales of the microparticle geometry, the aggregate scale geometry 
and the macroscopic scale of interest we have determined criteria for selecting which conditions 
are most applicable. We find that on small scales the error induced by a no-slip boundary 
condition on the aggregate surface is negligible. However, on larger scales these errors add 
up and the more accurate Beavers and Joseph condition or the Saffman approximation must 
be used. 

The applicability of equations (2.20) to curved aggregate surfaces is assumed. However, it 
is seen that for large micropore sizes this induces notable errors in the approximation. There 
is clearly scope for improving this approximation based on the curvature of the aggregate 
surface. The assumption of periodicity is applied to the internal structure of the aggregate, the 
aggregate surface and the aggregate scale packing. For real geometries obtained from X-ray 
CT, this assumption is not completely valid. Physically, the error induced by this assumption 
can be achieved by choosing a sufficiently large aggregate scale geometry which may include 
multiple aggregates. Validation of this assumption could be achieved by comparing several 
different regions of a single soil sample obtained from X-ray CT. This would not only validate the 
theory developed in this paper but also determine the required dimensions Ly and Lz required for 
accurate estimates of hydraulic conductivity to be obtained. 

The modelling in this paper was developed in the context of fluid flow in soil. However, 
it is applicable to a much wider set of situations in the study of porous media, for example 
petroleum reservoirs. This work highlights the importance of the role of the different scales in 
determining the macroscopic properties of fluid flow in porous media and how different aspects 
of the geometry contribute to these values. 
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